Geomarketing avec R - Session 2 - Boufarsi Adil

Introduction

Pour cette deuxième séance, nous allons apprendre à utiliser un planifateur d’itinéraire sur R afin de calculer des isochrones autour de chaque supermarché et ainsi pouvoir identifier les zones de chalandises.

1. Cartographier des isodistances

Une première information simple à obtenir mais très informative peut être obtenue en traçant des courbes d’isodistances autour de chaque supermarchés. On peut ainsi indentifier les zones non-deserviés.

• Représentez les courbes d’isodistances (st_buffer(), tm_borders())

télécharger le package en .zip en cliquant sur ce lien
https://cran.r-project.org/src/contrib/Archive/opentripplanner/opentripplanner_0.4.0.tar.gz

Et installer Java 8 sur votre ordinateur:
https://www.java.com/fr/download/

Show/Hide

setwd("C:/Users/boufarsi/Documents/thèse/Enseignement/Geomarketing/New/")

Avis : The working directory was changed to C:/Users/boufarsi/Documents/thèse/Enseignement/Geomarketing/New inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

pacman::p_load(dplyr, ggplot2, sf,readr,readxl,tmap,tmaptools)

#installer le package opentripplaner (en modifiant le chemin d'acces)
install.packages("C:/Users/boufarsi/Documents/thèse/Enseignement/Geomarketing/New/opentripplanner_0.4.0.tar.gz", repos = NULL, type = "source")

WARNING: Rtools is required to build R packages but is not currently installed. Please download and install the appropriate version of Rtools before proceeding:

https://cran.rstudio.com/bin/windows/Rtools/

Warning in install.packages :
  le package ‘opentripplanner’ est en cours d'utilisation et ne sera pas installé

library(opentripplanner)
iris_grenoble<-st_read("iris_grenoble.gpkg")

Reading layer `iris_grenoble' from data source `C:\Users\boufarsi\Documents\thèse\Enseignement\Geomarketing\New\iris_grenoble.gpkg' using driver `GPKG'
Simple feature collection with 70 features and 36 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 5.678578 ymin: 45.15415 xmax: 5.753078 ymax: 45.21432
Geodetic CRS:  WGS 84

#vérifiez bien que vous avez la derniere version du fichier, avec 22 lignes (22 points de vente)
sirene_grenoble<-st_read("sirene_grenoble.gpkg")

Reading layer `sirene_grenoble' from data source `C:\Users\boufarsi\Documents\thèse\Enseignement\Geomarketing\New\sirene_grenoble.gpkg' using driver `GPKG'
Simple feature collection with 21 features and 106 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 5.70661 ymin: 45.16637 xmax: 5.744077 ymax: 45.20106
Geodetic CRS:  WGS 84

st_crs(sirene_grenoble) #l'unité est le mètre

Coordinate Reference System:
  User input: WGS 84 
  wkt:
GEOGCRS["WGS 84",
    ENSEMBLE["World Geodetic System 1984 ensemble",
        MEMBER["World Geodetic System 1984 (Transit)"],
        MEMBER["World Geodetic System 1984 (G730)"],
        MEMBER["World Geodetic System 1984 (G873)"],
        MEMBER["World Geodetic System 1984 (G1150)"],
        MEMBER["World Geodetic System 1984 (G1674)"],
        MEMBER["World Geodetic System 1984 (G1762)"],
        MEMBER["World Geodetic System 1984 (G2139)"],
        ELLIPSOID["WGS 84",6378137,298.257223563,
            LENGTHUNIT["metre",1]],
        ENSEMBLEACCURACY[2.0]],
    PRIMEM["Greenwich",0,
        ANGLEUNIT["degree",0.0174532925199433]],
    CS[ellipsoidal,2],
        AXIS["geodetic latitude (Lat)",north,
            ORDER[1],
            ANGLEUNIT["degree",0.0174532925199433]],
        AXIS["geodetic longitude (Lon)",east,
            ORDER[2],
            ANGLEUNIT["degree",0.0174532925199433]],
    USAGE[
        SCOPE["Horizontal component of 3D system."],
        AREA["World."],
        BBOX[-90,-180,90,180]],
    ID["EPSG",4326]]

buffer<-st_buffer(sirene_grenoble,500) #200m isodistance

tmap_mode("view")
tm_shape(iris_grenoble) +tm_text(text = "iris_name",size=0.8)+
  tm_polygons("Pop_15_64", alpha=0.5,
              style="pretty", id="iris_name",
              title="Population (15-64 ans)")+
     tm_shape(buffer)+tm_borders("chartreuse3",lwd=3)+
  tm_shape(sirene_grenoble)+
  tm_dots("enseigne1et",legend.show = TRUE,title="Enseigne",id="enseigne1et",popup.vars="trancheeffe.2",alpha=0.3)+
  tm_layout(legend.outside = TRUE)

• Avec des cercles pleins plutôt que des courbes

Show/Hide

tm_shape(iris_grenoble) +tm_text(text = "iris_name",size=0.8)+
  tm_polygons("Pop_15_64", alpha=0.5,
              style="pretty", id="iris_name",
              title="Population (15-64 ans)")+
     tm_shape(buffer)+tm_polygons("enseigne1et",legend.show = FALSE,id="enseigne1et",alpha=0.5)+
   tm_shape(sirene_grenoble)+
  tm_dots("enseigne1et",legend.show = TRUE,title="Enseigne",id="enseigne1et",popup.vars="trancheeffe.2",alpha=0.8)+
  tm_layout(legend.outside = TRUE)

NA

2. Obtenir des isochrones

2.1 Mise en place du calculateur d’itinéraire

• Nous allons utiliser OpenTripPlanner, un logiciel codé en Java et utilisé notamment dans le calculateur d’itinéraire du TAG. Pour fonctionner, OTP a besoin de 3 choses:

Un dossier & un sous-dossier specifiques à créer dans votre R working directory > OTP > graphs > default

Des données sur le réseau de transport GTFS (General Transit Feed Specification) à déposer dans OTP > graphs > default

Un fond de carte openstreetmap à déposer dans OTP > graphs > default

• données GTFS

Show/Hide

knitr::include_graphics("C:/Users/boufarsi/Documents/thèse/Enseignement/Geomarketing/New/gtfs.png")

Show/Hide

Vous pouvez par exemple trouver les données GTFS sur le site M mobilité du TAG: data.mobilites-m.fr
Une fois téléchargé le fichier SEM-GTFS.zip renommez le en gtfs.zip et placez le dans OTP > graphs > default

• fond de carte OSM

Show/Hide

Il y a aussi beaucoup de sources, par exemple :
https://download.geofabrik.de/europe/france/rhone-alpes.html
Placez le fichier dans OTP > graphs > default

• Mise en place d’OTP

path_data <- "C:/Users/boufarsi/Documents/thèse/Enseignement/Geomarketing/New/OTP" #tell him the OTP directory 
path_otp <- otp_dl_jar() #java 

#Vous n'avez besoin de lancer cette ligne qu'une fois
#log <- otp_build_graph(otp = path_otp, # Build Graph
 #                      dir = path_data,memory=15000)
log1<-otp_setup(otp = path_otp, dir = path_data) # Start OTP
otpcon <- otp_connect()

Show/Hide

setwd("C:/Users/boufarsi/Documents/thèse/Enseignement/Geomarketing/New/")

Avis : The working directory was changed to C:/Users/boufarsi/Documents/thèse/Enseignement/Geomarketing/New inside a notebook chunk. The working directory will be reset when the chunk is finished running. Use the knitr root.dir option in the setup chunk to change the working directory for notebook chunks.

knitr::include_graphics("C:/Users/boufarsi/Documents/thèse/Enseignement/Geomarketing/New/OTP_capture.png")

2.2 Estimer les isochrones (otp_isochrone())

Show/Hide

iso<-otp_isochrone(otpcon,fromPlace = sirene_grenoble,fromID=sirene_grenoble$siret,mode = c("WALK","TRANSIT"),
            maxWalkDistance = 2000,date_time=as.POSIXct(strptime("2022-02-22 08:35", "%Y-%m-%d %H:%M")),
            cutoffSec = c(5, 7,9,11) * 60 ) # Cut offs in seconds

iso$minutes = iso$time / 60

#verification
colnames(iso)[3]<-"siret"
iso<-merge(iso,st_drop_geometry(sirene_grenoble[,c(3,43)]),by="siret",all.y=TRUE)%>%unique() #pour avoir les enseignes mais on a pas d'isochrone pour deux points de vente
sirene_grenoble[16,]$geom[[1]]<-st_point(c(5.728253270386639,45.17659752397455), dim = "XYZ") #ce point ne marchait pas pour 7min

#On recommence
iso<-otp_isochrone(otpcon,fromPlace = sirene_grenoble,fromID=sirene_grenoble$siret,mode = c("WALK","TRANSIT"),
                   maxWalkDistance = 1500,date_time=as.POSIXct(strptime("2022-02-22 08:35", "%Y-%m-%d %H:%M")),
                   cutoffSec = c(5, 7,9,11) * 60 ) # Cut offs in seconds
colnames(iso)[3]<-"siret"

iso<-merge(iso,st_drop_geometry(sirene_grenoble[,c(3,43)]),by="siret",all.y=TRUE)%>%unique() #pour avoir les enseignes mais on a pas d'isochrone pour deux points de vente
iso$minutes = iso$time / 60

3. Cartes finales et indices

3.1 Carte intéractive avec isochrones

Show/Hide

tmap_mode("view")
tm_shape(iris_grenoble) +#tm_text(text = "iris_name",size=0.8)+
  tm_polygons("Pop_15_64", alpha=0.5,
              style="pretty", id="iris_name",
              title="Population (15-64 ans)")+
   tm_shape(sirene_grenoble) + 
  tm_dots("enseigne1et",legend.show = TRUE,title="Enseigne",size=0.1,id="enseigne1et",popup.vars="trancheeffe.2")+
  tm_layout(legend.outside = TRUE)+
  tm_shape(iso) +  
  tm_fill("minutes",
          breaks = c(0, 5.01,  7.01,9.01,11.01), title="Isochrones (minutes)",
          style = "fixed",labels =c("0 to 5", "5 to 7", "7 to 9","9 to 11"),
          palette ="-BuPu",id="minutes",alpha = 0.3) +
  tm_borders()

3.2 Analyse de la concurrence

• Calculez le nombre de supermarchés accessibles par IRIS (st_intersects())

Show/Hide

iso420<-iso%>%filter(time==420) #on selectionne seulement l'isocrhone à 7min
#on compte (rowsums) le nombre de supermarchés présent par IRIS(st_intersect)
iris_grenoble$nb_supermarche_per_iris<-rowSums(ifelse(st_intersects(x=iris_grenoble,y=iso420,sparse=F),1,0)) 
tm_shape(iris_grenoble)+tm_polygons("nb_supermarche_per_iris",
              style="pretty", id="iris_name",labels =c("0", "1","2","3","4","5"),
              title="Nb de Supermarchés accessible en 7min")+tm_text(text = "iris_name",size=0.8)

• Bonus pour la 3ème session, on calcule le nombre de supermarchés de chaque enseigne accessibles par IRIS

Show/Hide

data_temp<-matrix(0,nrow=nrow(iris_grenoble),ncol=nrow(iso420))
for(j in 1:nrow(iso420)){
  for (i in 1:nrow(iris_grenoble)) {
    data_temp[i,j]<-as.character(st_intersects(iris_grenoble, iso420,sparse=F)[i,j])
    data_temp[i,j]<-ifelse(data_temp[i,j]=="TRUE",iso420$enseigne1et[j],0)
  }
} #cela nous donne les supermarchés accessibles dans chaque IRIS

data_temp2<-matrix(0,nrow=nrow(data_temp),ncol=length(unique(iso420$enseigne1et)))
for(j in 1:length(unique(iso420$enseigne1et))){
  data_temp2[,j]<-apply(data_temp, 1, function(x) length(which(x==unique(iso420$enseigne1et)[j])))
} #regroupe le dataframe différente avec le nombre de supermarché de chaque enseigne par IRIS

data_temp2<-as.data.frame(data_temp2)
colnames(data_temp2)<-unique(iso420$enseigne1et)
#pensez bien à changer le repertoire ici
write.csv(data_temp2,"C:/Users/boufarsi/Documents/thèse/Enseignement/Geomarketing/New/nb_pdv_iris.csv")

• Calculez la surface cumulée couverte par un supermarché pour chaque iris (st_intersection(), st_area())

Show/Hide

#Surface cumulée couverte par un supermarché par IRIS

intersect_pct <- st_intersection(iris_grenoble, iso420) %>% #calcul de l'intersection entre les iris et les iso
     mutate(intersect_area = st_area(.)) %>%   # création d'une variable de la surface en intersection
    dplyr::select(iris_name, intersect_area) %>%   # garder que les colonnes utiles
  group_by(iris_name)%>%mutate(intersect_area=sum(intersect_area))%>% #somme de la surface en intersection avec tous les supermarchés par iris
   st_drop_geometry()%>%unique() #On a pas besoin des informations géographiques

Avis : attribute variables are assumed to be spatially constant throughout all geometries

• Puis calculez la surface de chaque IRIS, le taux de couverture par iris et faîtes une carte

Show/Hide

#Surface total des IRIS
iris_grenoble <- mutate(iris_grenoble, iris_area = st_area(iris_grenoble)) #obtenir la surface de chaque iris

#Fusionner les deux 
iris_grenoble <- merge(iris_grenoble, intersect_pct, by = "iris_name", all.x = TRUE) #merge, which introduce some NAs when intersect_area=0

#Calcul du taux de couverture
iris_grenoble <- iris_grenoble %>% mutate(intersect_area=ifelse(is.na(intersect_area),0,intersect_area))%>% 
   mutate(couverture = as.numeric(intersect_area/iris_area)*100) #on obtient le taux de couverture

#Et enfin une carte
tm_shape(iris_grenoble)+tm_polygons(col="couverture", alpha=0.5,
              style="cont", id="iris_name",
              title="Couverture (%)")+tm_text(text = "iris_name",size=0.8)+
   tm_shape(sirene_grenoble) + 
  tm_dots("enseigne1et",legend.show = TRUE,title="Enseigne",size=0.1,id="enseigne1et",popup.vars="trancheeffe.2")+
  tm_layout(legend.outside = TRUE)

NA

• Enregistrer iris_grenoble pour la prochaine séance

st_write(iris_grenoble,"C:/Users/boufarsi/Documents/thèse/Enseignement/Geomarketing/New/iris_grenoble2.gpkg")

Layer iris_grenoble2 in dataset C:/Users/boufarsi/Documents/thèse/Enseignement/Geomarketing/New/iris_grenoble2.gpkg already exists:
use either append=TRUE to append to layer or append=FALSE to overwrite layer

Error: Dataset already exists.

Devoirs pour la prochaine séance

Préparez plusieurs cartes avec des isochrones différentes et d’autres variables au niveau des IRIS (carte 3.1) afin de proposer une localisation potentielle d’un nouveau supermarché à Grenoble. Proposez aussi une autre carte (3.2) construite à partir d’une isochrone différente (>7min). Ce travail sera évalué.